clear all; clc;close all;
%Updated: 14/08/215

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');
ratio=load('ratio.txt');
pg(:) =ratio(:,1);
GDP(:) =ratio(:,2);
pp(:) =ratio(:,3);

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');

na=rows(a);
nk=1;
ns=na/nk;
nb=rows(b);

alpha = 0.75;
pm     = 1;
gamma  = 0.43;
omegac = 0.3;
lambda = 0.184;
r      = 0.02;


atrans = ones(nb,1)*a';
lambda = 0.184;
r=0.02;

for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((i-1)*nb + j);
        
        p(i,j)  =p0((i-1)*nb + j);
        h(i,j)  =h0((i-1)*nb + j);
        ibp(i,j)=ibp0((i-1)*nb + j);
        ct(i,j) =ct0((i-1)*nb + j);
        ca(i,j) =ca0((i-1)*nb + j);
        gn(i,j) =gn0((i-1)*nb + j);       
        tau(i,j)=tau0((i-1)*nb + j);       
        w(i,j)  =w0((i-1)*nb + j);   
        def(i,j)=def0((i-1)*nb + j);   
        q(i,j)  =q0((i-1)*nb + j);   
        rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
        spr(i,j) =  (1+q(i,j).^(-1)-lambda)./(1+r)-1;

    end;
end;

cn =h.^alpha-gn;
const = (ct.^omegac).*(cn.^(1d0-omegac));
bp = b(ibp);

            
for i=1:na;
    vaut(i)  =vaut0(i);
    
    paut(i)  =paut0(i);
    haut(i)  =haut0(i);
    ctaut(i) =ctaut0(i);
    gnaut(i) =gnaut0(i);       
    tauaut(i)=tauaut0(i);       
    waut(i)  =waut0(i);   
end;

cnaut =haut.^alpha-gnaut;
constaut = (ctaut.^omegac).*(cnaut.^(1d0-omegac));

for i=1:na;
    temp1 = find(ct(i,:) > 0);
    temp2 = find(cn(i,:) > 0);
    temp3 = find(p(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    v(i,1:imin(i)-1)  = NaN; 
    p(i,1:imin(i)-1)  = NaN; 
    bp(i,1:imin(i)-1) = NaN; 
    ibp(i,1:imin(i)-1)= NaN; 
    ct(i,1:imin(i)-1) = NaN; 
    cn(i,1:imin(i)-1) = NaN; 
    const(i,1:imin(i)-1) = NaN;     
    gn(i,1:imin(i)-1) = NaN; 
    h(i,1:imin(i)-1)  = NaN; 
    ca(i,1:imin(i)-1)  = NaN; 
    tau(i,1:imin(i)-1)= NaN; 
	rec(i,1:imin(i)-1)= NaN;
%    def(i,1:imin(i)-1)= NaN; 
	q(i,1:imin(i)-1)  = NaN;

%     temp1 = find(ctaut(i,:) > 0);
%     temp2 = find(cnaut(i,:) > 0);
%     temp3 = find(paut(i,:) > 0);
% 
%     iminctaut(i) = temp1(1);
%     imincnaut(i) = temp2(1);
%     iminpaut(i)  = temp3(1);
% 
%     imin(i)  = max(iminctaut(i),max(imincnaut(i),iminpaut(i)));
% 
%     vaut(i,1:imin(i)-1)  = NaN; 
%     paut(i,1:imin(i)-1)  = NaN; 
%     ctaut(i,1:imin(i)-1) = NaN; 
%     cnaut(i,1:imin(i)-1) = NaN; 
%     constaut(i,1:imin(i)-1) = NaN;     
%     gnaut(i,1:imin(i)-1) = NaN; 
%     haut(i,1:imin(i)-1)  = NaN; 
%     maut(i,1:imin(i)-1)  = NaN; 
%     tauaut(i,1:imin(i)-1)= NaN; 
% 
end;

for i=1:ns;
for j=1:nk;        
    vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
    p2(i,j,:)  = p((i-1)*nk+j,:);  
    bp2(i,j,:) = bp((i-1)*nk+j,:);  
    ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
    ct2(i,j,:) = ct((i-1)*nk+j,:);  
    cn2(i,j,:) = cn((i-1)*nk+j,:);  
    const2(i,j,:) = const((i-1)*nk+j,:);      
    gn2(i,j,:) = gn((i-1)*nk+j,:);  
	ca2(i,j,:) = ca((i-1)*nk+j,:);  
    h2(i,j,:)  = h((i-1)*nk+j,:);   
    tau2(i,j,:)= tau((i-1)*nk+j,:);  
	rec2(i,j,:)= rec((i-1)*nk+j,:); 
    def2(i,j,:)= def((i-1)*nk+j,:);  
	q2(i,j,:)  = q((i-1)*nk+j,:); 
	ww2(i,j,:) = alpha*p2(i,j,:).*h2(i,j,:).^(alpha-1); 
    

    vaut2(i,j)  = vaut((i-1)*nk+j); 
    paut2(i,j)  = paut((i-1)*nk+j);  
    ctaut2(i,j) = ctaut((i-1)*nk+j);  
    cnaut2(i,j) = cnaut((i-1)*nk+j);  
    constaut2(i,j) = constaut((i-1)*nk+j);      
    gnaut2(i,j) = gnaut((i-1)*nk+j);  
    haut2(i,j)  = haut((i-1)*nk+j);  
    tauaut2(i,j)= tauaut((i-1)*nk+j);  
	wwaut2(i,j) = alpha*paut2(i,j).*haut2(i,j).^(alpha-1); 
end;
end;


lowa = round(2*ns/6);
meda = round(ns/2);
hgha = round(4*ns/6);

lambda=0.184;
r = 0.02;
b =-b/(0.184+0.02);  %change of sign

%default decision and prices
def21(:,:) = def2(:,1,:);
def22(:,:) = def2(:,1,:);

price21(:,:) = q2(:,1,:);
price22(:,:) = q2(:,1,:);

figure;
subplot(2,2,1);
contourf(b(1:nb),a(1:ns),def21(:,:),20);
xlabel('b');
ylabel('a');
title('default decision \kappa_{L}');
subplot(2,2,2);
contourf(b(1:nb),a(1:ns),def22(:,:),20);
xlabel('b');
ylabel('a');
title('default decision \kappa_{H}');
subplot(2,2,[3 4]);
plot(b(1:nb),price21(lowa,:),'-r',b(1:nb),price21(meda,:),':r',b(1:nb),price21(hgha,:),'--r',...
    b(1:nb),price22(lowa,:),'-b',b(1:nb),price22(meda,:),':b',b(1:nb),price22(hgha,:),'--b','LineWidth',1.5);
xlabel('b');


%%
%value functions: repayment vs autarky
vcon21(:,:) = vcon2(:,1,:);
vcon22(:,:) = vcon2(:,1,:);
vaut21(:,:) = vaut2(:,1);
vaut22(:,:) = vaut2(:,1);

figure;
plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
    vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
xlabel('b');title('value functions');
ylim([-16 -13]);legend('rep low y^T','rep high y^T','aut low y^T','aut high y^T')

%% Plot variables as function of debt for two yT levels
hr2  = h2;
cnr2 = cn2;
gnr2 = gn2;
ctr2 = ct2;
pr2  = p2;
taur2= tau2;
bpr2 = bp2;
wwr2 = ww2;
qr2  = q2;
spr2 = 100*((1+qr2.^(-1)-lambda)./(1+r)-1);
car2 = ca2;

%optimal policies and prices: repayment vs autarky
%NameArray = {'LineStyle','Color','LineWidth'};
%ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';

%NameArray0 = {'LineStyle','Color','LineWidth'};
%ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
NameArray = {'LineStyle','Color','LineWidth'};
ValueArray = {'--','--','-','-';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';

NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

[ilowa_aut]=find(def2(lowa,1,:)==0);
[ilowa_rep]=find(def2(lowa,1,:)==1);
hr2(lowa,1,ilowa_aut) = NaN;
ha2(lowa,1,1:nb) = haut2(lowa,1);
ha2(lowa,1,ilowa_rep) = NaN;

[ihgha_aut]=find(def2(hgha,1,:)==0);
[ihgha_rep]=find(def2(hgha,1,:)==1);
hr2(hgha,1,ihgha_aut) = NaN;
ha2(hgha,1,1:nb) = haut2(hgha,1);
ha2(hgha,1,ihgha_rep) = NaN;

pr2(lowa,1,ilowa_aut) = NaN;
pa2(lowa,1,1:nb) = paut2(lowa,1);
pa2(lowa,1,ilowa_rep) = NaN;

pr2(hgha,1,ihgha_aut) = NaN;
pa2(hgha,1,1:nb) = paut2(hgha,1);
pa2(hgha,1,ihgha_rep) = NaN;

gnr2(lowa,1,ilowa_aut) = NaN;
gna2(lowa,1,1:nb) = gnaut2(lowa,1);
gna2(lowa,1,ilowa_rep) = NaN;

gnr2(hgha,1,ihgha_aut) = NaN;
gna2(hgha,1,1:nb) = gnaut2(hgha,1);
gna2(hgha,1,ihgha_rep) = NaN;

bpr2(lowa,1,ilowa_aut) = NaN;
bpr2(hgha,1,ihgha_aut) = NaN;

cnr2(lowa,1,ilowa_aut) = NaN;
cna2(lowa,1,1:nb) = cnaut2(lowa,1);
cna2(lowa,1,ilowa_rep) = NaN;

cnr2(hgha,1,ihgha_aut) = NaN;
cna2(hgha,1,1:nb) = cnaut2(hgha,1);
cna2(hgha,1,ihgha_rep) = NaN;

ctr2(lowa,1,ilowa_aut) = NaN;
cta2(lowa,1,1:nb) = ctaut2(lowa,1);
cta2(lowa,1,ilowa_rep) = NaN;

ctr2(hgha,1,ihgha_aut) = NaN;
cta2(hgha,1,1:nb) = ctaut2(hgha,1);
cta2(hgha,1,ihgha_rep) = NaN;

wwr2(lowa,1,ilowa_aut) = NaN;
wwa2(lowa,1,1:nb) = wwaut2(lowa,1);
wwa2(lowa,1,ilowa_rep) = NaN;

wwr2(hgha,1,ihgha_aut) = NaN;
wwa2(hgha,1,1:nb) = wwaut2(hgha,1);
wwa2(hgha,1,ihgha_rep) = NaN;

taur2(lowa,1,ilowa_aut) = NaN;
taua2(lowa,1,1:nb) = tauaut2(lowa,1);
taua2(lowa,1,ilowa_rep) = NaN;

taur2(hgha,1,ihgha_aut) = NaN;
taua2(hgha,1,1:nb) = tauaut2(hgha,1);
taua2(hgha,1,ihgha_rep) = NaN;

car2(lowa,1,ilowa_aut) = NaN;
car2(hgha,1,ihgha_aut) = NaN;

qr2(lowa,1,ilowa_aut) = NaN;
qr2(hgha,1,ihgha_aut) = NaN;

spr2(lowa,1,ilowa_aut) = NaN;
spr2(hgha,1,ihgha_aut) = NaN;

endb = 2;

figure;
subplot(4,2,1);P = plot(b,100*squeeze(1-hr2(lowa,1,:)),b,100*squeeze(1-ha2(lowa,1,:)),b,100*squeeze(1-hr2(hgha,1,:)),b,100*squeeze(1-ha2(hgha,1,:)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:10:30);
axis(AX1, [min(b) endb -1 30]);
set(P,NameArray,ValueArray);
subplot(4,2,2);P = plot(b,squeeze(pr2(lowa,1,:)),b,squeeze(pr2(hgha,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(pa2(hgha,1,:)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low y^T','high y^T','Location','NorthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',3.5:0.5:6);
axis(AX1, [min(b) endb 3.5 6]);
set(P,NameArray,ValueArray2);
subplot(4,2,3);P = plot(b,squeeze(gnr2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gnr2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0.12:0.04:0.28);
axis(AX1, [min(b) endb 0.12 0.28]);
set(P,NameArray,ValueArray);
subplot(4,2,4);P = plot(b,-squeeze(bpr2(lowa,1,:))/(lambda+r),b,-squeeze(bpr2(hgha,1,:))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-1.6:0.8:1.6);
axis(AX1, [min(b) endb -1.6 1.6]);
set(P,NameArray0,ValueArray0);
subplot(4,2,5);P = plot(b,squeeze(cnr2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cnr2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.1:0.9);
axis(AX1, [min(b) endb 0.6 0.9]);
set(P,NameArray,ValueArray);
subplot(4,2,6);P = plot(b,squeeze(ctr2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ctr2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.2:1.4);
axis(AX1, [min(b) endb 0.6 1.4]);
set(P,NameArray,ValueArray);
subplot(4,2,7);P = plot(b,squeeze(wwr2(lowa,1,:)),b,squeeze(wwa2(lowa,1,:)),b,squeeze(wwr2(hgha,1,:)),b,squeeze(wwa2(hgha,1,:)));
title('w','Fontweight','normal');xlim([min(b) max(b)]);
AX1 = gca;box on;
set(AX1,'YTick',3:0.5:5);
axis(AX1, [min(b) endb 3 5]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(4,2,8);P = plot(b,squeeze(spr2(lowa,1,:)),b,squeeze(spr2(hgha,1,:)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:20:60);
axis(AX1, [min(b) endb -2 60]);
%version 2: with current account
% subplot(4,2,8);P = plot(b,-squeeze(car2(lowa,1,:)),b,-squeeze(car2(hgha,1,:)));
% %xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
% title('current account','Fontweight','normal');
% AX1 = gca;box on;
% set(AX1,'YTick',-0.25:0.15:0.35);
% axis(AX1, [min(b) endb -0.25 0.35]);
% set(P,NameArray0,ValueArray0);
dbstop;

%% plot time series
seriesp   = load('seriesp.txt');
seriesr   = load('seriesr.txt');
seriesh   = load('seriesh.txt');
seriesb   = load('seriesb.txt');
seriesct  = load('seriesct.txt');
seriesm   = 1; %load('seriesm.txt');
seriesgn  = load('seriesgn.txt');
seriestau = load('seriestau.txt');    
seriesa   = load('seriesa.txt');    
seriesbor = load('seriesbor.txt');    

n = length(seriesp);

seriescn= seriesh.^alpha-seriesgn;
seriesy = seriesa.*seriesm.^gamma+seriesp.*seriesh.^alpha; % -pm*seriesm;
seriesgny = seriesgn./seriesy;

time =1:n;

ValueArray1 = {'-';'blue';1.5}';
ValueArray2 = {'-','--';'red','blue';1.5,1.5}';
ValueArray3 = {'-','--','-';'red','blue','black';1.5,1.5,1.5}';

% No shaded areas for default episodes
% figure;
% subplot(3,2,1); P=plot(time,seriesa');title('a');
% set(P,NameArray,ValueArray1);
% subplot(3,2,2);P=plot(time,seriesp);title('p');
% set(P,NameArray,ValueArray1);
% subplot(3,2,3);P=plot(time,seriesct);
% %AX=legend('c_{T}','Location','NorthWest','Orientation','Horizontal');
% title('c_T');
% set(P,NameArray,ValueArray1);
% %LEG = findobj(AX,'type','text');
% %set(LEG,'FontSize',7);
% subplot(3,2,4); PPP=plot(time,seriescn,time,seriesh,time,seriesgn);
% axis([1 n -0.05 1.05]);
% AX=legend('h','c_{N}','g_{N}','Location','NorthWest','Orientation','Horizontal');title('h,c_{N} & g_{N}');
% set(PPP,NameArray,ValueArray3);
% LEG = findobj(AX,'type','text');
% set(LEG,'FontSize',7);
% subplot(3,2,5); P=plot(time,seriesr);title('spreads');
% set(P,NameArray,ValueArray1);axis([1 n -0.01 0.1]);
% subplot(3,2,6); P=plot(time,-seriesb);title('b');
% set(P,NameArray,ValueArray1);


seriesbor = 1-seriesbor;
seriesbor(199)=1;
seriesbor(200)=0;
seriesh(200) = 0;
seriesr(200) = 0.1;
ntime = 200;

figure;
subplot(3,2,1); 
shadedTimeSeries(time,seriesa',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('a');
AX1 = gca;box on;
set(AX1,'YTick',0.85:0.1:1.15);
axis(AX1, [time(1) time(ntime) 0.85 1.15]);
subplot(3,2,2);
shadedTimeSeries(time,seriesp',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('p');
AX1 = gca;box on;
set(AX1,'YTick',0:0.2:1);
axis(AX1, [time(1) time(ntime) 0.95*min(seriesp) 1.05*max(seriesp)]);
subplot(3,2,3);
shadedTimeSeries(time,seriesct',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('c_T');
AX1 = gca;box on;
set(AX1,'YTick',0.05:0.05:1.15);
axis(AX1, [time(1) time(ntime) 0.9*min(seriesct)  1.1*max(seriesct)]); 
subplot(3,2,4); 
shadedTimeSeries(time,seriesh',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('h,c_N & g_N');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
hold on;
H0 = line(time,seriesh','Color','b','LineWidth',2);
H1 = line(time,seriescn','Color','r','LineWidth',2,'LineStyle','--');
H2 = line(time,seriesgn','Color','g','LineWidth',2,'LineStyle',':');
hh = [H0, H1,H2];
legend(hh,'h','c_N','g_N','Location','SouthEast','Orientation','Horizontal');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
subplot(3,2,5); 
shadedTimeSeries(time,seriesr',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('spreads');
AX1 = gca;box on;
set(AX1,'YTick',0.0:0.025:0.1);
axis(AX1, [time(1) time(ntime) 0 0.1]);
subplot(3,2,6); P=plot(time,-seriesb);title('b');
shadedTimeSeries(time,-seriesb',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('b');
AX1 = gca;box on;
%set(AX1,'YTick',0.0:0.02:ceil(max(-1.1*seriesb/0.001))*0.001);
axis(AX1, [time(1) time(ntime) 0 max(-seriesb)]);

dbstop;